HEALPix logo

UXarray for Advanced HEALPix Analysis & Visualization

In this section, you’ll learn:

  • Using the uxarray package to perform advanced analysis operators over HEALPix data such as non-conservative zonal means, etc.

Prerequisites

Concepts

Importance

Notes

UXarray

Necessary

HEALPix overview

Necessary

Time to learn: 30 minutes

import uxarray as ux
import cartopy.crs as ccrs

Open data catalog

Tip

We assume, you have already gone over the previous section, UXarray for Basic HEALPix Statistics & Visualization. If not and if you need to learn about data catalogs in general and the data we will use throughout this notebook, we recommend to check that section first.:::

Let us open the online catalog from the WCRP’s Digital Earths Global Hackathon 2025 using intake and read the output of the ICON run ngc4008, which is stored in the HEALPix format:

import intake

# Final data catalog location (once hackathon website (https://digital-earths-global-hackathon.github.io/) updated)
#cat_url='https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml'
# Interim data catalog location
cat_url='https://raw.githubusercontent.com/digital-earths-global-hackathon/catalog/refs/heads/ncar/online/main.yaml'
cat = intake.open_catalog(cat_url)
model_run = cat.icon_ngc4008

We can look into a fine resolution dataset at zoome level = 10 in it as Xarray.Dataset:

ds = model_run(zoom=9, time="P1D").to_dask()
ds

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

<xarray.Dataset> Size: 232TB
Dimensions:                              (time: 10958, depth_half: 73,
                                          cell: 3145728, level_full: 90,
                                          crs: 1, depth_full: 72,
                                          soil_depth_water_level: 5,
                                          level_half: 91,
                                          soil_depth_energy_level: 5)
Coordinates:
  * crs                                  (crs) float32 4B nan
  * depth_full                           (depth_full) float32 288B 1.0 ... 5....
  * depth_half                           (depth_half) float32 292B 0.0 ... 5....
  * level_full                           (level_full) int32 360B 1 2 3 ... 89 90
  * level_half                           (level_half) int32 364B 1 2 3 ... 90 91
  * soil_depth_energy_level              (soil_depth_energy_level) float32 20B ...
  * soil_depth_water_level               (soil_depth_water_level) float32 20B ...
  * time                                 (time) datetime64[ns] 88kB 2020-01-0...
Dimensions without coordinates: cell
Data variables: (12/103)
    A_tracer_v_to                        (time, depth_half, cell) float32 10TB ...
    FrshFlux_IceSalt                     (time, cell) float32 138GB ...
    FrshFlux_TotalIce                    (time, cell) float32 138GB ...
    Qbot                                 (time, cell) float32 138GB ...
    Qtop                                 (time, cell) float32 138GB ...
    Wind_Speed_10m                       (time, cell) float32 138GB ...
    ...                                   ...
    vas                                  (time, cell) float32 138GB ...
    w                                    (time, depth_half, cell) float32 10TB ...
    wa_phy                               (time, level_half, cell) float32 13TB ...
    zg                                   (level_full, cell) float32 1GB ...
    zghalf                               (level_half, cell) float32 1GB ...
    zos                                  (time, cell) float32 138GB ...

Create UXarray Datasets from HEALPix

We can use from_healpix as follows to open a HEALPix grid from xarray.Dataset:

uxds = ux.UxDataset.from_healpix(ds)

uxds

<xarray.UxDataset> Size: 232TB
Dimensions:                              (time: 10958, depth_half: 73,
                                          n_face: 3145728, level_full: 90,
                                          crs: 1, depth_full: 72,
                                          soil_depth_water_level: 5,
                                          level_half: 91,
                                          soil_depth_energy_level: 5)
Coordinates:
  * crs                                  (crs) float32 4B nan
  * depth_full                           (depth_full) float32 288B 1.0 ... 5....
  * depth_half                           (depth_half) float32 292B 0.0 ... 5....
  * level_full                           (level_full) int32 360B 1 2 3 ... 89 90
  * level_half                           (level_half) int32 364B 1 2 3 ... 90 91
  * soil_depth_energy_level              (soil_depth_energy_level) float32 20B ...
  * soil_depth_water_level               (soil_depth_water_level) float32 20B ...
  * time                                 (time) datetime64[ns] 88kB 2020-01-0...
Dimensions without coordinates: n_face
Data variables: (12/103)
    A_tracer_v_to                        (time, depth_half, n_face) float32 10TB ...
    FrshFlux_IceSalt                     (time, n_face) float32 138GB ...
    FrshFlux_TotalIce                    (time, n_face) float32 138GB ...
    Qbot                                 (time, n_face) float32 138GB ...
    Qtop                                 (time, n_face) float32 138GB ...
    Wind_Speed_10m                       (time, n_face) float32 138GB ...
    ...                                   ...
    vas                                  (time, n_face) float32 138GB ...
    w                                    (time, depth_half, n_face) float32 10TB ...
    wa_phy                               (time, level_half, n_face) float32 13TB ...
    zg                                   (level_full, n_face) float32 1GB ...
    zghalf                               (level_half, n_face) float32 1GB ...
    zos                                  (time, n_face) float32 138GB ...
uxds.uxgrid.face_node_connectivity

<xarray.DataArray 'face_node_connectivity' (n_face: 3145728, n_max_face_nodes: 4)> Size: 101MB
array([[ 778523, 1700461, 2064568, 1448407],
       [ 261422, 2488188, 1700461,  778523],
       [1700461, 1778639, 2171953, 2064568],
       ...,
       [1147106,  293699,  791866,  409789],
       [ 791866, 2110736, 1740665, 2548812],
       [ 293699, 2286345, 2110736,  791866]], shape=(3145728, 4))
Dimensions without coordinates: n_face, n_max_face_nodes
Attributes:
    cf_role:      face_node_connectivity
    long name:    Maps every face to its corner nodes.
    start_index:  0
    _FillValue:   -9223372036854775808

Data variable of interest

Then let us pick a variable from the dataset, which will give us an uxarray.UxDataArray:

uxda = uxds['tas']
uxda

<xarray.UxDataArray 'tas' (time: 10958, n_face: 3145728)> Size: 138GB
[34470887424 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
Dimensions without coordinates: n_face
Attributes:
    cell_methods:   time: mean cell: mean
    component:      atmo
    grid_mapping:   crs
    long_name:      temperature in 2m
    standard_name:  air_temperature
    units:          K
    vgrid:          height_2m

Global mean and plot

Computing the global air temperature mean (at the first timestep) and also having a quick plot of it would be a good idea to have as references to compare the upcoming analyses & visualizations to them:

%%time
print("Global air temperature average on ", uxda.time[0].values, ": ", uxda.isel(time=0).mean().values, " K")

Global air temperature average on  2020-01-02T00:00:00.000000000 :  286.3096  K
CPU times: user 490 ms, sys: 279 ms, total: 769 ms
Wall time: 1.61 s

%%time

projection = ccrs.Robinson()

uxda.isel(time=0).plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global temperature",
    width=700,
)

CPU times: user 12.4 s, sys: 588 ms, total: 13 s
Wall time: 13.6 s

WARNING:param.GeoOverlayPlot00477: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

Cross-sections

We can look at constant latitude/longitude cross-sections of an uxarray.UxDataArray:

boulder_lat = 40.0190


# With fine resolutions like zoom level of 9, it is visually hard to observe the cross-sections, 
# so we will use a zoom level of 4 for a better visualization
uxda_coarse = ux.UxDataset.from_healpix(model_run(zoom=4, time="P1D").to_dask())['tas']
uxda_coarse.uxgrid.face_node_connectivity

uxda_lat = uxda_coarse.cross_section.constant_latitude(boulder_lat)
uxda_lat

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

<xarray.UxDataArray 'tas' (time: 10958, n_face: 128)> Size: 6MB
[1402624 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
Dimensions without coordinates: n_face
Attributes:
    cell_methods:   time: mean cell: mean
    component:      atmo
    grid_mapping:   crs
    long_name:      temperature in 2m
    standard_name:  air_temperature
    units:          K
    vgrid:          height_2m
import geoviews.feature as gf

uxda_lat.isel(time=0).plot(
    rasterize=False,
    projection=projection,
    global_extent=True,
    cmap="inferno",
    clim=(220, 310),
    features=["coastline"],
    title=f"Global temperature cross-section at {boulder_lat} degrees latitude",
    width=700,
) * gf.grid(projection=projection)

/home/runner/miniconda3/envs/healpix-cookbook-dev/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_graticules_30.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

WARNING:param.GeoOverlayPlot00713: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

Let’s also look at the mean of the cross-section:

print(f"Mean at {boulder_lat} degrees lat (Boulder, CO, USA): {uxda_lat.mean().values} K")

Mean at 40.019 degrees lat (Boulder, CO, USA): 286.6632080078125 K

Latitude interval

uxda_lat_interval = uxda_coarse.cross_section.constant_latitude_interval([boulder_lat-15, boulder_lat+15])

uxda_lat_interval.isel(time=0).plot(
    rasterize=False,
    projection=projection,
    global_extent=True,
    cmap="inferno",
    clim=(220, 310),
    features=["coastline"],
    title=f"Global temperature cross-section at the latitude interval [{boulder_lat-5},{boulder_lat+5}] degrees",
    width=700,
) * gf.grid(projection=projection)

WARNING:param.GeoOverlayPlot00937: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

print(f"Mean at the latitude interval of [{boulder_lat-5},{boulder_lat+5}] degrees (-/+15 degrees Boulder, CO, USA): {uxda_lat_interval.mean().values} K")

Mean at the latitude interval of [35.019,45.019] degrees (-/+15 degrees Boulder, CO, USA): 286.2819519042969 K

Non-conservative Zonal mean

Calculating the zonal mean is easy by providing the latitude range between -90 and 90 degrees with a step size in degrees:

zonal_mean_tas = uxda.isel(time=0).zonal_mean(lat=(-90, 90, 5))

(
    uxda.isel(time=0).plot(
        cmap="inferno",
        # periodic_elements="split",
        height=300,
        width=600,
        colorbar=False,
        ylim=(-90, 90),
    )
    + zonal_mean_tas.plot.line(
        x="tas_zonal_mean",
        y="latitudes",
        height=300,
        width=180,
        ylabel="",
        ylim=(-90, 90),
        xlim=(220, 310),
        # xticks=[220, 250, 280, 310],
        yticks=[-90, -45, 0, 45, 90],
        grid=True,
    )
).opts(title="Temperature and its Zonal means at every 5 degrees latitude")

Remapping

Now, we will be looking into a remapping case. The data set we are using in this section has the zoom=10 available but not for all the variables, e.g. tas has all NaN values at that zoom level. Let us try to remap our data here that is sampled at zoom=9 to a Grid at zoom=10.

Let’s start with creating the destination uxarray.Grid:

%%time
uxgrid_zoom10 = ux.Grid.from_healpix(zoom=10, pixels_only=False)

CPU times: user 20.5 s, sys: 4.58 s, total: 25.1 s
Wall time: 10.4 s

%%time 
uxda_remapped = uxda.isel(time=0).remap.inverse_distance_weighted(
    uxgrid_zoom10, k=3, remap_to="face centers"
)

CPU times: user 8min 51s, sys: 752 ms, total: 8min 51s
Wall time: 8min 52s

%%time

uxda_remapped.plot(
    projection=projection,
    cmap="inferno",
    features=["borders", "coastline"],
    title="Global temperature - remapped to zoom=10",
    width=700,
)

CPU times: user 35.6 s, sys: 1.81 s, total: 37.5 s
Wall time: 37.8 s

WARNING:param.GeoOverlayPlot01467: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
On this page